﻿using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using System.Threading.Tasks;
using System.Threading;
using System.Collections;

namespace Dwarf.Managed
{
    /// <summary>
    /// Presentation of Hidden Markov model.
    /// </summary>
    class HMM
    {
        /// <summary>
        /// Number of states;  Q = {1,2,...,N} .
        /// </summary>
        private int numberOfStates;

        /// <summary>
        /// Number of observation symbols; V={1,2,...,M}.
        /// </summary>
        private int numberOfObservationSymbols;

        /// <summary>
        /// A[1..N][1..N]. a[i][j] is the transition prob
        /// of going from state i at time t to state j
        /// at time t+1.
        /// </summary>
        public double[][] stateTransitionMatrix;

        /// <summary>
        /// B[1..N][1..M] is the  probability 
        /// of observing symbol j in state i .
        /// </summary>
        public double[][] observationProbabilityMatrix;
        /// <summary>
        /// [1..numberOfStates]  Is the initial state distribution.
        /// </summary>
        public double[] initialStates;

        public HMM(int n, int m)
        {
            numberOfStates = n;
            numberOfObservationSymbols = m;
            //creation of  state transition and observation matrices
            stateTransitionMatrix = new double[numberOfStates][];
            observationProbabilityMatrix = new double[numberOfStates][];
            for (int i = 0; i < numberOfStates; i++)
            {
                stateTransitionMatrix[i] = new double[numberOfStates];
                observationProbabilityMatrix[i] = new double[numberOfObservationSymbols];
            }

            initialStates = new double[numberOfStates];
        }

        public int getNumberOfStates()
        {
            return numberOfStates;
        }
        public int getNumberOfObservation()
        {
            return numberOfObservationSymbols;
        }
    }

    /// <summary>
    /// Presentation of Viterbi model.
    /// </summary>
    class ViterbiModel
    {
        /// <summary>
        /// Length of observation the sequence.
        /// </summary>
        private int lengthOfObservationSequence;
        /// <summary>
        /// Maximization of probability.
        /// </summary>
        private double probability;

        /// <summary>
        /// Is the best score(highest probability) along a single path,at time t.
        /// </summary>
        public double[][] delta;

        /// <summary>
        /// Track of the argument which maximized.
        /// </summary>
        public int[][] psi;

        /// <summary>
        /// Observation sequence [ O ].
        /// </summary>
        public int[] observationSequence;

        /// <summary>
        /// Sequence of state [ Q ].
        /// </summary>
        public int[] sequenceOfState;

        public ViterbiModel(int t, HMM hmm)
        {
            lengthOfObservationSequence = t;
            //creation of delta and psi work elements
            delta = new double[lengthOfObservationSequence][];
            psi = new int[lengthOfObservationSequence][];
            int numberOfState = hmm.getNumberOfStates();
            for (int i = 0; i < lengthOfObservationSequence; i++)
            {
                delta[i] = new double[numberOfState];
                psi[i] = new int[numberOfState];
            }

            observationSequence = new int[lengthOfObservationSequence];
            sequenceOfState = new int[lengthOfObservationSequence];
        }

        public int getLengthOfObservationSequence()
        {
            return lengthOfObservationSequence;
        }
        public double getProbability()
        {
            return probability;
        }
        public void setProbability(double p)
        {
            probability = p;
        }
    }

    /// <summary>
    /// Solve the Graph model problem.
    /// </summary>
    /// <remarks> 
    /// For the demonstration purposes of the project,
    /// we have chosen Viterbi algorithm for Hidden Markov Model
    /// </remarks>
    class Solver
    {
        /// <summary>
        /// Instance of Hidden Markov model.
        /// </summary>
        public HMM hmm;

        /// <summary>
        /// Instance of Viterbi model.
        /// </summary>
        public ViterbiModel vit;

        public Solver()
        {
        }

        /// <summary>
        /// Based method.
        /// </summary>
        public void Solve()
        {

            /* 1. Initialization  delta and psi */
            int numberOfStates = hmm.getNumberOfStates();
            Parallel.For(0, numberOfStates, delegate(int i)
            {
                vit.delta[0][i] = hmm.initialStates[i] * (hmm.observationProbabilityMatrix[i][vit.observationSequence[0]]);
                vit.psi[0][i] = 0;
            });

            /* 2. Compute */
            int lengthOfObsSeq = vit.getLengthOfObservationSequence();
            for (int t = 1; t < lengthOfObsSeq; t++)
            {
                Parallel.For(0, numberOfStates, delegate(int j)
                {
                    double maxval = 0.0;            //temporary storage for max value for delta array
                    int maxvalind = 0;              //temporary storage for max value for psi array
                    double val = 0.0;
                    //finding the maximum
                    for (int i = 0; i < numberOfStates; i++)
                    {
                        val = vit.delta[t - 1][i] * (hmm.stateTransitionMatrix[i][j]);
                        if (val > maxval)
                        {
                            maxval = val;
                            maxvalind = i;
                        }
                    }
                    vit.delta[t][j] = maxval * (hmm.observationProbabilityMatrix[j][vit.observationSequence[t]]);
                    vit.psi[t][j] = maxvalind;
                });
            }

            /* 3. Termination */

            vit.setProbability(0.0);

            vit.sequenceOfState[lengthOfObsSeq - 1] = 0;
            Parallel.For(0, numberOfStates, delegate(int i)
            {
                if (vit.delta[lengthOfObsSeq - 1][i] > vit.getProbability())
                {
                    vit.setProbability(vit.delta[lengthOfObsSeq - 1][i]);
                    vit.sequenceOfState[lengthOfObsSeq - 1] = i;
                }
            });

            /* 4. Path (state sequence) backtracking */
            for (int t = lengthOfObsSeq - 2; t >= 0; t--)
            {
                vit.sequenceOfState[t] = vit.psi[t + 1][vit.sequenceOfState[t + 1]];
            }
        }

    }
}